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Abstract 

A C++ class was written for the calculation of frequentist confidence intervals using 
the profile likelihood method. Seven combinations of Binomial, Gaussian, Poissonian 
and Binomial uncertainties are implemented. The package provides routines for the 
calculation of upper and lower limits, sensitivity and related properties. It also 
supports hypothesis tests which take uncertainties into account. It can be used in 
compiled C++ code, in Python or interactively via the ROOT analysis framework. 
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• Title of Program: TRolke version 2.0 

• Program available from: CPC Program Library, ... 

• Licensing provisions: MIT License 

• Computer for which the program is designed: Unix, GNU/Linux, Mac 

• Operating Systems under which the program has been tested: Linux 2.6 (Sci- 
entific Linux 4 and 5, Ubuntu 8.10) , Darwin 9.0 (Mac-OS X 10.5.8) 

• Programming Language used: ISO C++ 

• Memory required to execute with typical data: ~ 20 MB, 

• No. of bytes in distributed program, including initialization file, etc.. 1 MB 

• Distribution Format: tar file 

• Keywords: confidence interval calculation, systematic uncertainties, profile 
likelihood 

• Nature of the Physical Problem: The problem is to calculate a frequentist 
confidence interval on the parameter of a Poisson process with statistical or 
systematic uncertainties in signal efficiency or background. 

• Method of solution: Profile likelihood method, Analytical 

• Typical Running Time: < 10~ 4 seconds per extracted limit. 
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1 Introduction and scope 



Routines were written for the calculation of frequentist confidence intervals 
using the profile likelihood method. The package provides routines for the 
calculation of upper and lower limits, average limits (sensitivity) and related 
properties, taking uncertainties in background estimate and signal efficiency 
into account. The implementation considers seven different statistical models 
with different combinations of Binomial, Gaussian, Poissonian or no uncertain- 
ties. For example in the Gaussian background case, our package derives upper 
and lower limits on the signal strength for a Poisson process with Gaussian 
background expectation b ± 5b. It is also possible to construct hypothesis tests 
which take uncertainties into account. The statistical problems are treated 
using the Profile Likelihood method. 

The package provides a C + + class with accompanying examples. It can be 
used in compiled code, interactively via the ROOT [1J analysis framework, 
and from Python. This is TRolke version 2.0. It adds to version 1 (imple- 
mented in Fortran and in C++): hypothesis tests, a reworked user interface, 
documentation, examples and python support. 

This paper is organized as follows. First, the profile likelihood methods is 
summarized, section [2] ; second, it is shown how our routines can be used 
for optimization of statistical discovery or limit setting power, section [3j The 
means for specification of the statistical model, and in general the class inter- 
face are described in section HI 
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2 The profile likelihood method 

Frequentist limits are constructed from data such that when repeated with 
new data the limits cover the fixed but unknown parameter value it with a 
frequency which converges to the requested probability, the confidence level 
1 — a. Limit calculation methods are often based on the inversion of an hypoth- 
esis test, as described in e.g. [2] [3] [3] • an d we follow the same scheme. Classical 
hypothesis tests investigate the validity of a default hypothesis, the null hy- 
pothesis Hq; that an examined sample of data is compatible with background 
and we call the complementary hypothesis T-l-y a discovery. The profile likeli- 
hood method is based on the likelihood ratio tests statistic now described. For 
some observable X, let us assume a probability density function /(Xj|7r, b) 
depending on k parameters 7r = {71*1, . . . , 7Tfc} of interest to the researcher (such 
as the strengths of different signal sources), and I additional nuisance param- 
eters b = {b\, . . . , bi} (such as the strength of different background sources). 
For a set of n independent observations X = {X\, . . . , X n } the likelihood is 



where the denominator is the likelihood maximized over the whole {tt, b} 
space, while the nominator is maximized over the more restrictive null hy- 
pothesis space {77 = 7To, b}. The likelihood ratio A is also known as the profile 
likelihood and is a stochastic function explicitly depending on the data (and 
the null hypothesis) but not the nuisance parameters. 



n 



L(7r,b\X) = l[f(X l \7r,b). 



i=i 



The likelihood ratio test statistic is defined as 




SUp{L(7T, b\X); 77 = 7r , b} 

sup{L(7r, b\X); 7T, b} 
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In general the inversion of a test to find the confidence region requires scan- 
ning over all possible signals, as described for example in [2]. Our routines 
instead make use of a very powerful result from mathematical statistics, that 
under some general conditions the distribution of —2 log A converges to a chi- 
square distribution with k degrees of freedom. Although these conditions are 
not satisfied in the problem considered here it has been shown that its per- 
formance is surprisingly good, especially when, as here, nuisance parameters 
are included. The statistical performance of the Profile likelihood method is 
studied in Ref.JS]. 



3 Analysis optimization for optimal limits or discovery power 

In this section we describe how our routines are used for optimization of anal- 
ysis cuts, with the figure of merit being either stringent limits (in case the 
signal is expected to be weak), or probability for discovery (if the signal is 
expected to be strong). 

3. 1 Analysis optimization for stringent limits in case of vanshing signal 

When a signal is expected to be weak enough so that significant discovery 
is unlikely, it is relevant to optimize the analysis for optimal limit setting 
power. This can be done by assuming no signal and minimizing the so-called 
sensitivity. For example with a 90% confidence level (that is, a =10%), let 
us denote a calculated upper limit sgo- The sensitivity of the experiment is 
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defined as the average upper limit in case of vanishing signal; 



s 9oO) = P ( X i b ) S 90(x, b), (1) 

2 = 



where P(x, b) is the Poisson probability of observing x events for background 
expectation b, in absence of signal. For finding the optimal analysis cut we can 
assume without loss of generality that the background and signal expectations 
are monotonically decreasing functions of a cut c: s(c) = /i s e s (c), and 6(c). The 
constant /i s is the assumed normalisation of the signal at some arbitrary "no 
cut" level so that all uncertainties in the signal rate expectation are attributed 
to the signal detection efficiency e s . 

As an example, let's consider an energy dependent spectrum of particles 
probed by a particle detector. For the physical test spectrum 

dE =° test ^^' (2) 
the expected number of observed signal events is 

nest = "test ' T J dfi j ^±a(E)dE. (3) 

The cross section a determines the detection efficiency which is now a function 
of the energy E, and T is the exposure time. For the observation of x events 
the model rejection factor £ (x) is defined as 

£(x) = s 90 ( y x ,b)/s test . (4) 

The upper limit can be written in terms of the test signal 
d^ (E) _ dgtestCg) 

~ltE~ ~ m ' dE ' (5j 
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and the average limit on the signal strength, set by repeated independent 
experiments in case of vanishing signal is 

d^9o(£) _ T d$ test (E) 

dE dE ' {> 

where £ = sgo(&)/ s test * s called the model rejection potential. 

Our package provides s 90 through the method GetSensitivity(L>& s^, D&i sjj), 
and the upper limit sgo (and the lower limits) through Getl_imits(Z}& s-^, D&i 
sjj), where D indicates a double precision value. 

3.2 Hypothesis testing with uncertainties 

In order to reject "Hq with significance a, the number of observed events 
Xq must be equal to or higher than a critical number Xc(b), where b is the 
background expectation. The significance a is the probability of observing x c 
or more events from a stochastic background with mean b assuming vanishing 
signal. 

The part of sample space rejecting "Hq is called the critical region, while its 
complement is called the acceptance region. In the constructed test the crit- 
ical region is completely defined by x c . If the background expectation b was 
completely known, we could find xq by solving 

oo 

P(n > x c \b) = £ P(n\b) < a, (7) 

n=x c 

where P(n\b) is the Poisson distribution, but in general the background ex- 
pectation is unknown and so we find the critical value by inverting the profile 
likelihood method. Remembering that confidence regions are constructed such 
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that the true but unknown signal strength 5 is outside the confidence region 
with probability a for any fixed 5 we assume the hypothesis Hq which means 
5 = 0. The critical region is therefore defined as the subset of values x which 
gives rise to limits not covering 5 = 0. That is, Hq is rejected for observations 
that lead to lower limits sl larger than zero. The limits are monotonic in x, so 
the hypothesis test is completely characterised by a critical number x c , and 
written x > xq. This critical number algorithm is implemented as the method 
GetCriticalNumber(int& n c ). 

3.3 Analysis optimization for signal discovery 

Assuming a specific signal strength s = 5, it is relevant to consider the prob- 
ability of making a discovery. This is given by the power of the hypothesis 
test, F/3 = 1 — f3. A signal hypothesis H- s ^ is said to be at the visibility 
threshold if it leads to a discovery with a pre-specified probability Fp, for ex- 
ample 50%. Discovery is claimed when xo > xc, so in order to minimise the 
visibility threshold, signal is added to the (background) expectation until the 
probability for x > x<z is at least Fp. 

For the case of vanishing uncertainties, the visibility threshold can be directly 
calculated [B] from the Poisson distribution by finding the smallest signal s^p 
fulfilling 

P(n>a; c p|6 + s th p)>i5> (8) 

or equivalently P{n < x c p|6 + s^p) < 0, where the critical value x c p is 
that found using equation [71 The quantity s^p is the visibility threshold for 
the signal expectation in case of vanishing uncertainties. The construction is 
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shown in figure [TJ Uncertainties are accounted for through the critical num- 
0.30i 1 




Fig. 1. For a predefined /3, the visibility threshold is the smallest signal that is 
discovered with at least probability Fp = 1 — /3 at significance a. In this example, 
a = 1%, (3 = 50%, b = 3.5, x c = 8, s th = 4.17. 

ber xc(ct, b, A&) as function of significance and expectation number. A method 
similar to this has previously been described by PunziJT]. As in equation [HJ sig- 
nal is added to the (background) expectation until the probability for rejection 
of Hq is at least Fp. This means 

x c (Fp, b + s, A b+S ) > x c {a, b, A b ) (9) 

where A& and Ab+ S represent the total uncertainties of background, and back- 
ground plus signal respectively. Equation [9] is solved numerically by finding 
the smallest allowed signal expectation s and the solution is called s^. Since 
the tested hypothesis "Hq assumes exactly S = 0, we do not include any uncer- 
tainty in the signal efficiency, while here the background estimate is assumed 
Gaussian. The described procedure for finding the critical number in the pres- 
ence of uncertainties is thus a function on the form xc(cv, b, A&), where A^ is 
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the background uncertainty. 

Assuming that the uncertainties of signal efficiency and the background esti- 
mate are sufficiently uncorrelated and Gaussian (or exact), equation [Hlbecomes 

x c (Fp, 6(c) + 8 (c), ^/A b (c) 2 + A s (c) 2 ) > x c (a, 6(c), A 6 (c)) . (10) 

For the observation of x events the model rejection factor £ (x) is defined as 
f(x) = s 9 o(x ,6)/s test , (11) 



where stest i s > as m section |373| the expectation number of signal events for an 
assumed test signal. The optimal cut c and the corresponding critical number 
Xq is found by minimising the signal strength as function of the cut c. 
The visibility threshold for the expected number of observed signal events for 
a fixed cut c is 

s th = ^ th e s (c)- (12) 

The physical threshold signal strength is found in terms of the test spectrum 
(in analogy with equation EJ) by = a^ es ^ • Sth/ S test> or equivalently 

d$ th(£) _ „ dgtegtXg) (u] 
dE ~ V ' dE ' 

where rj = s^/stest * s the model detection potential. Minimizing 77 optimizes 
the analysis such that the signal strength required for detection (with at least 
probability Fp — 1 — f$) is minimized. Our code provides the critical number 
and the through GetCriticalNumber(int& n c ) and 
bool TRolke2::Getl_eastDetectableSignal( J D& s th , D 
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4 Class interface and use 

The library allows seven combinations of efficiency and background rate mod- 
els, each presented here. Once the model and its parameters are specified, 
the user can obtain limits, critical numbers and so on as explained in the 
subsequent sections. 



4-1 Model Specification methods 

4-1.1 SetGaussBkgGaussEfF(x, bm, em, sde, sdb) 
Background: Gaussian, Efficiency: Gaussian 

This model implements the case of Gaussian background with expectation 
bm and standard deviation sdb and Gaussian efficiency with expectation 
em and standard deviation sde. The integer x is the number of observed 
events. 

4-1.2 SetGaussBkgKnownEff(x, bm, sdb, e) 
Background: Gaussian, Efficiency: Known 

This model implements the case of Gaussian background with expectation 
bm and standard deviation sdb and known efficiency e. The integer x is the 
number of observed events. 
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4-1.3 SetKnownBkgGaussEff(x, em, sde, b) 



Background: Known, Efficiency: Gaussian 

This model implements the case of Gaussian efficiency with expectation em 
and standard deviation sde and known background b. The integer x is the 
number of observed events. 



4-1-4 SetKnownBkgBinomEff(x, z, b, m) 



Background: Known, Efficiency: Binomial 

This model implements the case of known background expectation b and 
Binomial signal efficiency. The integer z is the number of observed events 
(in the signal region) out of the m evaluated signal (Monte Carlo) events. 
The integer x is the number of observed events. 
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4-1.5 SetPoissonBkgKnownEfF(x, y, r, e) 



Background: Poisson, Efficiency: Known 

The background is either measured simultaneously with signal, from side- 
bands, or with separate background Monte Carlo. The real value r is the 
size of the background region in terms of the size of the background regions. 
It can be used in two ways - Either it's the ratio between the size of the 
background and the signal regions in case background is observed (from 
sidebands), or in case background is determined from simulations; the ra- 
tio between simulated and observed exposure time. The background in the 
signal region is estimated from r and the integer y, the number of observed 
events in background region. The integer x is the number of observed events; 
as always in the signal region. 



4-1.6 SetPoissonBkgBinomEfF(x, y, z, r, m) 



Background: Poisson, Efficiency: Binomial 

This model implements the case of Binomial signal efficiency and Poissonian 
background estimate. For an explanation of Binomial efficiencies, please 
refer to section 14.1.41 and for Poissonian backgrounds to section 14.1.51 The 
integer x is the number of observed events. 
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4-1.7 SetPoissonBkgGaussEfF(x, y, em, sde, r) 



Background: Poisson, Efficiency: Gaussian 

This model implements the case of Gaussian signal efficiency and Poissonian 
background estimate. For an explanation of Binomial efficiency, please refer 
to I4.1.6[ and for Poissonian backgrounds to section 14.1.41 The integer x is 
the number of observed events. 

4-2 Configuration methods and constructor 

The confidence level (CL) is set either at object construction via an optional 
argument or with either of the SetCL or SetCLSigmas methods. 

Two options are offered to deal with cases where the maximum likelihood 
estimate (MLE) is not in the physical region. Bounding is controled with 
the SetBounding method. The "bounded likelihood" option corresponds to the 
"bounds for the physical region" option in MINUIT/MINOSjgJjS]. Unbounded 
likelihood allows the maximum likelihood estimate to be in the unphysical 
region. It has better coverage [5] and is used by default. 

4-3 Limit calculation methods 

The calculation of limits for the model and parameters as specified, is per- 
formed with the any of the following methods; 
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4.3.1 bool GetLimits(D& s L , D& stf 



This method calculates and returns the upper and lower limits for the pre- 
specified model, confidence level and model parameters. 

4.3.2 bool GetSensitivity(D& s L ,D& sjj) 

This method returns the average upper and average lower limits assuming 
vanishing signal. The summation is a Poisson sum over the background 
expectation. This can be used for cut optimization as described in sectioin 

4.3.3 bool GetLimitsQuantile(D& s L , D& sjj, int& out_x, D q = 0.5) 

This method returns the upper and lower limits for the outcome correspond- 
ing to a given quantile q assuming vanishing signal and a simple Poisson 
summation using the background expectation. As a default, the quantile 
value 0.5 is used, corresponding to median limits. The quantile and median 
method has the advantage over the sensitivity that it is independent of the 
signal parameter metric. The quantile x value is returned as out_x 

4.3.4 bool GetLimitsML(D& s L , D& sjj, int& out_x) 

This method provides the upper and lower limits for the most likely outcome 
(out_x), assuming vanishing signal. 
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4-4 Hypothesis test methods 



These two methods are used for hypothesis testing as described in section 13.31 

4-4-1 bool GetCriticalNumber(int& n c ) 

Get the smallest number of observed events x, corresponding to rejection of 
the null hypothesis. 

44.2 bool TRolke2::GetLeastDetectableSignal(D& s fh , D (3) 

Get the smallest signal strength leading to rejection of the null hypothesis 
with probability /3 as described in section 13.21 Currently Gaussian as well 
as vanishing uncertainties are supported. 

4-5 Availability and prerequisites 

The latest versions of the code, its documentation and examples are freely 
available [10J. The class makes use of a number of ROOT [TJ routines for stan- 
dard mathematical functions, the interactive interface and bindings which 
makes it easy to use our methods in Python. Examples of all functionality of 
the C++ class are included in our code and demonstrate its use with Python, 
as interactive C++, and as a compiled example program. 



17 



References 



[1] Rene Brun and Fons Rademakers, ROOT - An Object Oriented Data Analysis 
Framework, Proceedings AIHENP'96 Workshop, Lausanne, Sep. 1996, Nucl. 
Inst. & Meth. in Phys. Res. A 389 (1997) 81-86. 
http:/ /root.cern.ch/ 

[2] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57 (1998) 3873 



|arXiv:physics/9711021|. 



[3] A. Stuart and J. K. Ord: Kendall's Advanced Theory of Statistics, Vol. 2, 
Classical Inference and Relationship, Oxford University Press, New York (1991). 

[4] J. Neyman, Phil. Trans. Royal Soc. London A, 333, (1937). 

[5] W.Rolke, A. Lopez, J. Conrad, Nucl.Instrum.Meth.A551:493-503,2005 
|arXiv:hep-ph/0403059] . 



[6] G. C. Hill, Phys. Rev. D 67 (2003) 118101 arXiv:physics/0302057 



[7] G Punzi, Sensitivity of searches for new signals and its optimization, 
PHYSTAT2003 Sep8-ll 

[8] F. James, M. Roos, MINUIT, a System for Function Minimization and Analysis 
of the Parameter Errors and Correlations, Comput. Phys. Commun .10 (1975) 
343-367 

[9] F. James, Interpretation of the Shape of the Likelihood Function around its 
Minimum, Comput. Phys. Commun, Volume 20, Issue 1, (1980) 29-35 

[10] J. Lundberg, J. Conrad, W. Rolke, A. Lopez, TRolke 2.0, 



http: / / cpc.cs.qub.ac.uk / summaries / AEFT_vl_0.html 



18 



